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We present a study of site and bond percolation on periodic lattices with 3 nearest neighbors per 
site. We have considered 3 lattices, with different symmetries, different underlying Bravais lattices, 
and different degrees of longer-range connections. As expected, we find that the site and bond 
percolation thresholds in all of the 3-connected lattices studied here are significantly higher than in 
diamond. Interestingly, thresholds for different lattices are similar to within a few percent, despite 
the differences between the lattices at scales beyond nearest and next-nearest neighbors. 



I. INTRODUCTION 

Percolation is one of the simplest phase transitions 
known: sites on a lattice are occupied at random un- 
til there is a path that can be traversed from one end 
of the lattice to the other, traveling only neighbor-to- 
neighbor on occupied sites [lj. The set of sites along this 
path or connected to this path (via neighbor-to-neighbor 
steps along occupied sites) is called the spanning clus- 
ter. In the limit of a large system size (linear dimension 
L ^> 1), the probability of forming a spanning cluster 
via random occupation of sites goes to zero below a criti- 
cal occupation probability per site p c and unity above p c . 
Percolation models have been used to study a great many 
phenomena, from transport in porous media to biological 
systems to gelation [2j. 

A well-known trend in percolation on periodic 3D lat- 
tices is that p c increases as the coordination number z 
(number of nearest neighbors per site) decreases [3| . (See 
Table U near the end for examples and more references.) 
This trend arises from the fact that, when there are 
fewer neighbors per site, there are fewer ways to navi- 
gate around an obstacle. Consequently, more sites have 
to be occupied to ensure a path spanning from one end 
of the system to the other. To our knowledge, the lowest 
coordination number studied for percolation on simple 
3D lattices is z = 4, the most prominent example being 
the diamond lattice with p c « 0.43 for site percolation [4] 
and « 0.388 for bond percolation |5|. 

However, z = 4 is not the lowest possible coordina- 
tion number. The lowest possible non- trivial coordina- 
tion number is z = 3[6|, and this coordination number 
has been realized in interesting families of 3D lattices, 
namely the (n, 3) families of lattices studied by A. F. 
Wells [3] . The 3 refers to the number of nearest neigh- 
bors (bonds) per site, and n is the smallest number of 
steps that one would have to take along the lattice sites 
to return to the same point. In lattice families with mul- 
tiple members, (n, 3) is followed by a letter, e.g. (10, 3)-a, 
(10,3)-b, etc. A total of 30 lattices have been identified 
in these families, with n values of 7, 8, 9, 10, and 12. 



The (n, 3) lattices offer a chance to study the interplay 
between nearest neighbors, higher- level connectivity (via 
n), and other aspects of lattice geometry (e.g. the differ- 
ences between, say, (10, 3)-a and (10, 3)-b) for the lowest 
non-trivial z value. We are only aware of one study of 
percolation on lattices in this family, namely the (10, 3)-a 
lattice [8], and that study focused on invasion percolation 
and trapping, rather than the standard site and bond 
percolation problems. 

Recently, the (10,3)-a lattice has attracted additional 
attention because of its unusually high symmetry, pos- 
sessing a property known as "strong isotropy." [9| Only 
one other 3D lattice (diamond), shares this property with 
(10,3)-a. The (10,3)-a structure is observed in a num- 
ber of interesting systems, e.g. block copolymers fioj. 
and it has been proposed as a possible structure for a 
metastable phase of carbon [11|. Motivated by this re- 
cent attention to (10, 3)-a and related lattices, as well as 
the general question of how high the percolation thresh- 
old can be for a 3D lattice with z = 3, we have studied 3 
lattices in this family: (10,3)-a, its close relative (10,3)- 
b, and (8,3)-a. These particular lattices span multiple 
n values, and are especially easy to study because they 
can be realized with bonds of uniform length and 120° 
bond angles, making them easy to construct with ball 
and stick models. 

In what follows, we first discuss the geometries of 
our three representative lattices in their simplest forms: 
equal bond lengths and 120° bond angles. We illustrate 
deformations of the lattices that preserve the topology 
(connections between nearest neighbors) but enable a 
mapping onto a cubic lattice, for convenience in enu- 
merating sites in a calculation. We then summarize key 
features of the Newman-Ziff algorithm[12j used to deter- 
mine the site and bond percolation thresholds of the lat- 
tices, and present the computed percolation thresholds. 
Finally, we compare our results with other lattices. 



II. THE LATTICES UNDER STUDY 



A. The (10, 3)-a lattice 
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a 4-atom basis. If we work in a coordinate system where 
the sites of the bcc lattice are at the corners of a cube, 
one site is at the origin (0,0,0), and the edges are of 
unit length and parallel to the standard x, and z axes, 
then the 3 other atoms in the basis are at (0, — 
(— |, |, 0), and Q, 0, — |). The atoms correspond to sites 
through 3 in the figure (i.e. the atoms at the origin 
and its 3 nearest neighbors), they sit in the (111) plane, 
and the bonds to them are separated by 120° angles. 

In percolation theory, the geometry of the bonds is less 
important than the presence of the bonds. A deformation 
of the lattice and bonds that preserves the connections 
between sites will not change the percolation threshold of 
the lattice. In the Newman-Ziff algorithm that we used, 
all of the relevant information on the lattice is stored in 
the "nn" array, which contains the labels of the nearest 
neighbors bonded to each site in the lattice. For conve- 
nience in enumerating sites, we have deformed the lattice 
to fit onto a simple cubic grid, while preserving the bond- 
ing structure of the lattice. Our numbering system for 
sites, and the deformations used to fit them onto a cubic 
grid, are shown in Fig. [U^b) and Fig. H^c). 




FIG. 1. (a) The (10,3)-a lattice, with sites numbered, (b) and 
(c) The deformations used to map the lattice onto a cubic grid 
for computational convenience. 



B. The (10,3)-b lattice 

The (10,3)-b lattice is unusual, because even if we spec- 
ify unit bond lengths and 120 degree bond angles we still 
have an unconstrained structural degree of freedom: Sup- 
pose that we define the positions of all of the lattice sites 
and also the basis. We can then deform the lattice in a 
manner that uniformly changes the spacing between lat- 
tice planes, while also displacing the atoms in the lattice 
planes, without changing any bond lengths or bond an- 
gles. We have chosen to describe this lattice in a form 
that maximizes the symmetry, as this form lends itself 
to easier visualization with ball-and-stick models. In this 
form, the lattice is body-centered tetragonal, with lattice 

vectors (>/3,0,0), (0,^,0), and (^,^,3^). 

The 4-atom basis is more complicated. Instead of a 
central atom with 3 neighbors, the basis is a chain of 

4 atoms, located at (0,0,0), (o, (o, ^,§),and 

( — "^'^)- These correspond to sites through 3 in 
the figure. As before, we also deform this lattice to rep- 
resent it on a simple cubic grid. The numbering system 
and deformation steps are illustrated in Fig.[2fb) and (c). 




FIG. 2. (a) The (10,3)-b lattice, with sites numbered, (b) and 
(c) The deformations used to map the lattice onto a cubic grid 
for computational convenience 
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C. The (8,3)-a lattice 



We realized the (8,3)-a lattice, shown in Fig. [3ja), in 
a structure that has he xago nal symmetry, w ith lattice 
vectors (-5/2, V3/6, 2^/2/3), (5/2,y/3/6,2y/2/S), and 
(0,4^/3/3, y/2/3). The x axis is the axis of hexagonal 
symmetry. 

The (8, 3)-a lattice has a six-atom basis, 
with unit bond lengths and atoms located at 
(0,0,0), (0,-V3/3,-V273), (1/2,-5^5/6,-^/273), 
(3/2,-5^/3/6,-^/273), (2,-v^3/3,- v / 273), and 
(2,0,0). These correspond to sites through 5 in the 
figure. The steps used to deform the (8,3)-a lattice so 
that it fits onto a cubic grid are illustrated in Fig. [3jb) 
and (c). 



(a) 




FIG. 3. (a) The (8,3)-a lattice, with sites numbered, (b) and 
(c) The deformations used to map the lattice onto a cubic 
grid for computational convenience 



III. SIMULATIONS 
A. Description of the algorithm 

We used the Newman- Ziff algorithm for identifying the 
percolation threshold of finite-sized clusters[12]. In brief, 
the algorithm works by occupying sites (or bonds) on a 
lattice of N one-at-a-time in a random order. Relation- 
ships between sites are defined in the array nn, which 
indicates the site numbers for the nearest neighbors of 
each site. At each step, after a site is occupied we check 
whether the occupation of the nth site produces a clus- 
ter that wraps around the entire system. Each cluster 
is assigned a "pointer" to a root (or parent) site, corre- 
sponding to the first site in the cluster, and also a vector 
that points in the direction of the parent site. As clusters 
grow and merge, pointers and vectors are updated to the 
root of the largest cluster involved in the merger. Wrap- 
ping is detected when a newly-occupied site (1) joins two 
portions of the same cluster and (2) the vectors going 
from the joined portions of the cluster to the root differ 
by something other than the displacement vector between 
the sites. We consider wrapping between parallel faces 
of the computational volume {i.e. along the x, and z 
directions) as well as diagonal wrapping. (See the paper 
by Newman and Ziff [12] for more details.) The iteration 
ends when wrapping is detected. Another lattice is then 
generated, and the process is repeated, until Nl lattices 
(typically 10 3 in our work) have been generated. Bond 
percolation is handled in a completely analogous manner, 
substituting bonds for sites. Two bonds are considered 
neighbors if they touch the same site. 

This process generates a plot of wrapping probabil- 
ity Rl vs. occupation fraction <\> = n/N, where L refers 
to the linear dimension of the lattice. In order to de- 
termine as a function of occupation probability (the 
usual quantity of interest in percolation theory) , it is nec- 
essary to convolve Rl with a binomial distribution: 

R l(p) = E ( n ) P n ^-P) N - n Rdn) (1) 

n=0 ^ ' 

Eq. (pQ) amounts to a weighted sum over all possible re- 
alizations of an TV-site lattice with site occupation prob- 
ability p. The number of realizations that wrap with n 
occupied sites enters the sum via Rl{ti). Different occu- 
pation numbers n will have different likelihoods of being 
realized at the same occupation probability p, i. e. differ- 
ent degeneracies, and this enters the sum via the binomial 
distributions. Occupation numbers that are not close to 
N • p are unlikely, and hence get low weight from the 
binomial distribution. 

To implement the algorithm we closely followed the 
code example given in the paper by Newman and Ziff, 
implementing it in Python 2.7 on a quad-core Intel pro- 
cessor in Linux (Ubuntu 12.04). We checked our code 
by determining the site and bond percolation thresholds 
of the 2D square, 2D honeycomb, and 3D simple cubic 
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lattices. Adapting the code to the lattices under study 
only required modification to the function that identifies 
the nearest neighbors of each site, as well as the vectors 
from each site to its neighbors. To validate the output for 
the lattices under study, we generated lattices with small 
numbers of sites (« 30) and had the code output the list 
of occupied sites when wrapping occurred. We verified 
by hand that (1) there was a cluster that wrapped the 
lattice, and (2) the removal of the most recently occupied 
site would cause the cluster to not wrap. 

Given a plot of the wrapping probability Rl(p) for dif- 
ferent lattice sizes L, it is possible to obtain an estimate 
of the percolation threshold p c by looking for the point 
where the plot crosses over from low wrapping probabil- 
ity to high (e.g. the steepest point on the plot). An 
example for the (10,3)-a lattice is shown in Fig. EJ For 
sufficiently large system sizes, this estimate of p c can be 
quite accurate. More efficient approaches can obtain high 
precision and accuracy by comparing Rl(p) plots for sev- 
eral different system sizes L. One common way of com- 
paring plots at different sizes is to use the scaling relation 
\p c (L) — p c \ oc L~ u ~ v ', where p c (L) is obtained from the 
cross-over of an Rl(p) plot and uo and v are universal 
exponents that depend only on dimension [l3|. 

However, a simpler approach, one that enables a very 
intuitive estimate of p c and the associated uncertainty, 
is to make plots with Rl(p) (for fixed p) on the vertical 
axis and L on the horizontal axis. For p > p c , Rl(p) is 
an increasing function of L, and for p < p c Rl(p) is a 
decreasing function of L. This follows directly from the 
fact that Rl(p) becomes more "step-like" as L increases. 
At p = p c , Rl(p) is independent of Z/[l4|. The value 
of p c can thus be estimated simply by looking at the 
different plots and identifying the flattest one. When the 
plot of Rl oscillated somewhat as a function of L (due to 
randomness in the simulations), we identified the value of 
p for which Rl oscillates without a pronounced upward 
or downward trend. 



B. The uncertainty in the percolation threshold 

There is an easy way to determine the uncertainty in 
p c based only on the simulation output, and without 
any assumptions about critical exponents. In a plot of 
Rl(p), the percolation transition is identified by look- 
ing for a value of p at which the wrapping probability 
is approximately stationary as the system size changes. 
However, the wrapping probabilities are determined from 
finite samples in Monte Carlo simulations, and (as de- 
scribed below) this gives rise to statistical fluctuations in 
Rl for all p values. A fluctuation SRl shifts the curve 
along the horizontal axis by an amount 5Rl/ ((IRl/ dp), 
where dRL/dp is the slope of the Rl vs. p curve. 

In order to get the uncertainty in Rl(p), let us first 
consider ^ asa function of occupation number n rather 
than occupation probability p. We use generate Nl lat- 
tices with n occupied sites, and use the Newman- Ziff al- 
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FIG. 4. (a) Wrapping probability Rl vs. occupation probabil- 
ity p for sites on the (10,3)-a lattice. The curves for different 
system sizes L cross around p c ~ 0.5714. (b) Zoom of previ- 
ous plot. 



gorithm to determine the number Nl * Rl(^) that wrap. 
If we were to do this process repeatedly, we would find 
that the number wrapping obeys a binomial distribution 
with mean NlRl(^) and variance NlRl(^)(^ — Rl(^)), 
and so the fraction Rl(ti) that wrap has mean Rl(ti) 
and variance Rl(ti)(1 — Rl(ti))/Nl. 

When we go from Rl(ti) to Rl(p), the convolu- 
tion with a binomial distribution means that Rl(p) is 
a weighted sum of different i?z,(n) values. However, 
Rl(ti + 1) is not statistically independent of Rl(ti), since 
it is generated by occupying one more site (or bond) 
in each of the lattices used to determine Rl(u). Near 
Pc Rl(k 7^ Npc) is approximately a linear function of 
n - Np c : 

R L (n) « R L (Np c ) + c * (n - Np c ) + 77(71) (2) 

where c is an unknown constant and 77 is a noise term. 

This expression for Rl(ti) is convolved with a binomial 
probability distribution that is approximately Gaussian 
for large N (and hence even about p c ). The normaliza- 
tion of the probability distribution means that the convo- 
lution with the first term gives Rl(Np c ). The symmetry 
of the peak of the distribution means that the convolu- 
tions with the second and third terms vanish. We thus 
conclude that Rl(p = Pc) ~ Rl(^ = Np c ), and so we 

will use J Rl ^ Rl ^> for the uncertainty in Rl(Pc)- This, 
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in turn gives the uncertainty in p c as: 

\ Rl(1-R l ) ldR L 
SPc = i N L I ~dp~ (3) 

We have neglected the variance of the weighted sum of 
the noise terms r](n) in Eq. (j2j), but that variance is small 
because 77(71) exhibits correlations for nearby n values: 
The fact that Rl (n) is an increasing function of n means 
that downward fluctuations of Rl{^) are limited in mag- 
nitude by the fluctuations of Rl{tl — 1). The fact that 
Rl < 1 reduces the probability of successive upward fluc- 
tuations. 

In Eq. (|3j), we use the slope of the steepest Rl vs. p 
curve. If the steepest curve were perfectly vertical, fluc- 
tuations of the other curves would be completely irrele- 
vant, and the point of intersection would remain on that 
curve at the value of p where it rises. Consequently, the 
finite slope of the steepest curve is the limiting factor in 
our determination of p c . 



IV. RESULTS 

Figures [5] through [7] show Rl vs. L (for different p 
values), for site and bond percolation on the lattices un- 
der study. The number of unit cells in a realization of 
each lattice is L 3 , so that the number of lattice sites 
is 4L 3 (for the (10,3)-a and (10,3)-b lattices) or 6L 3 
(for (8,3)-a). The number of bonds is 6L 3 (for (10,3)- 
a and (10,3)-b) or 9L 3 (for (8,3)-a). In all plots, site 
occupation probability p was incremented in steps of 
1 /number of sites or bonds for smallest L value. 

TABLE I. Site and bond percolation thresholds for important 
3-dimensional lattices with different coordination numbers z. 
Uncertainties are given in parentheses, and refer to the last 
one or two digits, depending on the number of digits in the 
uncertainty. Bibliographic references are in brackets []. 



lattice 


z 


p c (site) 


p c (bond) 


(8,3)-a 


3 


0.577962(33) 


0.555700(22) 


(10,3)-a 


3 


0.571404(40) 


0.551060(37) 


(10,3)-b 


3 


0.565442(40) 


0.546694(33) 


diamond 


4 


0.4301(4)fl5] 


0.3893(2) [15] 


simple cubic 


6 


0.3116080(4) [16] 


0.2488126(5) [17] 


bcc 


8 


0.2459615(10) [16] 


0.1802875(10) [17] 


fee 


12 


0.1992365(10)fl7] 


0.1201635(10)fl6] 


hep 


12 


0.1992555(10) [18] 


0.1201640(10) [18] 



a For a more exhaustive table of percolation thresholds in 
different systems, Prof. Robert Ziff regularly updates the 
Wikipedia page "Percolation Threshold." A snapshot from the 
date of this writing has been archived at 
|http : //www . webcitat ion . org/6CTpDz4BX 
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FIG. 5. Wrapping probability Rl vs. L for (a) site and (b) 
bond percolation on (10,3)-a lattices, for different occupation 
probabilities p. (Selected p values shown to right of plots.) 
The percolation threshold p c is the value of p that gives the 
flattest overall trend. The uncertainty is determined via Eq. 

V. DISCUSSION 

The percolation thresholds identified for the 3- 
coordinated lattices considered here are higher than typi- 
cal values for 3-dimensional lattices that have been stud- 
ied previously. This point is illustrated in Table HJ which 
shows percolation thresholds for a variety of common 
3-dimensional lattices, organized by their coordination 
number z. It is clear from the table that p c increases 
as z decreases. This makes intuitive sense: With lower 
coordination numbers, it is easier to destroy a spanning 
cluster by removing a few sites or bonds at key points, 
while in lattices with higher coordination numbers there 
are more paths that can be navigated to circumvent a 
missing site or bond. 

The information shown in Table [I] can be compared 
with an analytical expression |3j for the approximate de- 
pendence of p c on spatial dimension d and coordination 
number z: 

p c =p ((d-l)(z-l))- a d b (4) 

where a = 0.6160 for sites and 0.9346 bonds, b = for 
sites and b = a for bonds, and po = 1.2868 for sites and 
0.7541 for bonds. This comparison is done in Fig. [8] For 
ease of presentation we only show the (10,3)-a lattice on 
these plots, but the percolation thresholds of the other 
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FIG. 6. Wrapping probability Rl vs. L for (a) site and (b) 
bond percolation on (10,3)-b lattices, for different occupation 
probabilities p. (Selected p values shown to right of plots.) 
The percolation threshold p c is the value of p that gives the 
flattest overall trend. The uncertainty is determined via Eq. 



3-connected lattices differ by only about 1%, and would 
overlap the (10,3)-a lattice if included. Crucially, the 
percolation thresholds of the 3-connected lattices are very 
close to the theoretical plot. 

We have not explored the entire family of 3-connected 
nets studied by Wells. Thus, it is an open question as 
to whether there are simple 3-dimensional periodic lat- 
tices with higher percolation thresholds. We qualify the 
previous sentence with the word "simple" because it is 
intuitively obvious that one could often trivially increase 
the percolation threshold of a crystalline structure by 
inserting chains of 2-connected sites between sites with 
higher coordination numbers. The interesting question 
is whether one can construct crystals with higher p c 
without making use of 2-connected sites between higher- 
coordinated cites. 



In conclusion, we have used Monte Carlo simulations to 
determine the percolation thresholds of several interest- 
ing lattices that have not, to the best of our knowledge, 
been studied previously. We find that these lattices have 
substantially higher percolation thresholds than other 3- 
dimensional lattices, due to their low coordination num- 
bers. The results for both bond and percolation thresh- 
olds are very close to theoretical predictions that rely on 
coordination number and dimensionality. 
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FIG. 7. Wrapping probability Rl vs. L for (a) site and (b) 
bond percolation on (8,3)-a lattices, for different occupation 
probabilities p. (Selected p values shown to right of plots.) 
Although there is some oscillation, the percolation threshold 
p c is the value of p that gives the flattest overall trend. The 
uncertainty is determined via Eq. ©. 
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